! 壳单元主程序(未完成),real64位
program pdyna
    use iso_fortran_env , only: real64
    use reader_mod
    use writer_mod
    use shell_mod
    use sconstants_mod
    use ls_mod,only: do_shell_static_solve
    implicit none
    integer                  :: i
    integer                  :: iostat, nnum, dbc_num, fbc_num, amps_num, fileunit
    character(len=250)       :: nodefp,elemfp,dbcfp,fbcfp,ampfp,configfp,inpfile
    real(real64),allocatable :: nodes(:,:),dbc_val(:),fbc_val(:), amps_data(:,:)
    integer,allocatable      :: el_conn(:,:),dbc_pos(:,:),fbc_pos(:,:),fbc_amp(:)
    ! 求解设置变量 
    character(len=2)         :: solve_type
    real(real64)             :: inc_time,step_time,E,t,nu,rho
    real(real64)             :: alpha,beta  ! 瑞利阻尼系数
    integer                  :: nev,gdof
    logical                  :: opened
    ! 结果数组
    real(real64),allocatable :: nodal_u(:,:),nodal_f(:,:)
    
    inc_time    = 0.0_real64
    step_time   = 0.0_real64
    E           = 0.0_real64
    nu          = 0.0_real64
    t           = 0.0_real64
    rho         = 0.0_real64
    alpha       = 0.0_real64
    beta        = 0.0_real64
    nev         = 0
    ! ==========================================================
    !                       读取输入文件的数据
    ! ==========================================================
    inquire(unit=100,opened=opened)
    if (.not. opened) then
        call get_command_argument(number=1,value=inpfile) ! 获取命令行参数
        open(unit=100,status='old',file=inpfile,action='read',iostat=iostat)
        if (iostat /=0) then
            print *, 'open commandline args file failed'
            stop
        end if
        read(100,'(A)') nodefp      ! 读取节点坐标文件名
        read(100,'(A)') elemfp      ! 读取单元连接文件名
        read(100,'(A)') dbcfp       ! 读取dbc文件名
        read(100,'(A)') fbcfp       ! 读取fbc文件名
        read(100,'(A)') ampfp       ! 读取amp文件名
        read(100,'(A)') configfp    ! 读取配置文件名    
    end if
    close(100) ! 关闭文件

    print *, '==> Start reading input data !'
    ! 读取 node.txt 数据
    call load_nodes(nodefp,nodes,nnum,iostat)
    if (iostat > 0) write(*,*) "load node failed, iostat=", iostat
    ! 读取 element.txt 数据
    call load_elements(elemfp,4,el_conn,iostat)
    if (iostat > 0) write(*,*) "load element failed, iostat=", iostat
    
    ! 读取 dbc.txt 数据
    call load_dbcs(dbcfp,dbc_pos,dbc_val,dbc_num,iostat)
    if (iostat > 0) write(*,*) "load dbc failed, iostat=", iostat
    
    ! 读取 fbc.txt 数据
    call load_fbcs(fbcfp,fbc_pos,fbc_amp,fbc_val,fbc_num,iostat)
    if (iostat > 0) write(*,*) "load fbc failed, iostat=", iostat
    
    ! 读取 amp.txt 数据
    call load_amps(ampfp,amps_data,amps_num,iostat)
    if (iostat > 0) write(*,*) "load amps failed, iostat=", iostat

    ! 读取 config.txt 数据
    open(newunit=fileunit,file=configfp,status='old',action='read',iostat=iostat)
    if (iostat /=0) then
        print *, 'open config file failed'
        stop
    end if
    read(fileunit,'(A2)') solve_type                             ! 读取分析步类型

    ! -----------------------------------------------------------
    !   show infomations
    ! -----------------------------------------------------------
    print *,   '------------------Short description------------------'
    write(*,*) "node number:", nnum
    write(*,*) 'element number:', size(el_conn,1)
    write(*,*) 'dispalcement boundary condition number:', dbc_num
    write(*,*) 'force boundary condition number:', fbc_num
    write(*,*) 'amplititude curves number:', amps_num
    write(*,*) 'solve type:',solve_type
    print *,   '------------------------------------------------------'

    select case(solve_type)
        case(linear_static_step)
            read(fileunit,*) inc_time, step_time                ! 读取增量步时间
            read(fileunit,*) E, nu, t                           ! 读取材料参数,壳单元厚度
            close(fileunit)                                     ! 关闭配置文件

            print *, "==>Start linear static analysis"
            allocate(nodal_u(nnum,node_dof_num))                ! 结果数组分配内存
            allocate(nodal_f(nnum,node_dof_num))                ! 结果数组分配内存
            ! 调用线性分析子程序
            gdof = nnum*node_dof_num
            call do_shell_static_solve(nodes,el_conn,e,nu,t,dbc_pos,dbc_val,fbc_pos,fbc_val,gdof,nodal_u,nodal_f)
            call ls_writer(1,nodes,el_conn,nodal_u,nodal_f)     ! 结果输出
            ! 写入tecplot结果文件
            call get_command_argument(number=2, value=inpfile)  ! 获取结果文件名
            open(newunit=fileunit, file=inpfile, status='replace',action='write',iostat=iostat)
            if (iostat/=0) print *, 'open commandline args file failed'
            write(fileunit,'(A)') 'TITLE = "Linear Static Analysis Result"'
            write(fileunit,'(A)') 'VARIABLES = "X","Y","Z","Ux","Uy","Uz","RFx","RFy","RFz"'
            write(fileunit,'(A,I8,A,I10,A)') 'ZONE T="SHELL_MESH",N=',nnum,',E=',size(el_conn,1),',ZONETYPE=FEQuad,DATAPACKING=POINT'
            do i=1,nnum                                         ! 写入节点和节点数据
                write(fileunit,'(3F28.10,3F30.16,3F30.16)') nodes(i,1),nodes(i,2),nodes(i,3),&
                nodal_u(i,1),nodal_u(i,2),nodal_u(i,3),nodal_f(i,1),nodal_f(i,2),nodal_f(i,3)
            end do
            do i=1,size(el_conn,1)                              ! 写入单元连接定义
                write(fileunit,'(4I16)') el_conn(i,:)
            end do
            close(fileunit)
            deallocate(nodal_u,nodal_f)
            print *, "==>End linear static step"

        case(natural_frequency_step)
            write(*,*) "processing natural frequency step"
        case(transient_dynamics_step)
            write(*,*) "processing transient dynamics step"
    end select
end program pdyna
